DVF Script 4

1 Les données DVF à une échelle locale : les exemples de l’agglomération de La Rochelle

Ce script décrit plusieurs fonctions de manipulation de la dimension spatiale des données DVF. Il propose aussi plusieurs formes de représentation cartographique de ces données à une échelle de l’EPCI (La Rochelle agglomération).

1.1 Chargement des packages nécessaires

library(tidyverse) # Manipulation de données
library(sf) # Manipulation de données spatiales
library(mapsf) # Cartographie thématique
library(cartography) # Cartographie thématique
library(linemap) # cartes en pics
library(spdep) # pour calculer l'auto-corrélation spatiale
library(geoR) # pour calculer le semi-variogramme empirique
library(spatstat) # pour produire des surfaces lissées
library(terra) # pour le traitement de données matricielles (raster)
library(cluster) # pour la CAH
library(happign) #pour récupérer les référentiels de l'IGN
library(maptiles) #pour ajouter des fonds aux représentattions carto
library(stars) #pour convertir des objets ppp
library(raster) #pour découper un raster par l'emprise d'un vecteur

1.2 Import du jeu de données brut

DVF <- read.csv("data/dvf_rochelle.csv", encoding = "UTF-8")

1.2.1 Création d’un tableau avec les indicateurs immobiliers de LR

RecapLR <- DVF %>% group_by(type) %>% 
  summarise(tot = n(), prixmed = round(median(prix)), prixmoy = round(mean(prix)), surfmed = round(median(surface)), surfmoy = round(mean(surface)), prixm2med = round(median(prix_m2)), prixm2moy = round(mean(prix_m2)))

RecapLR <- RecapLR %>% mutate(part = (tot/sum(tot)*100))
print(RecapLR)

## # A tibble: 2 × 9
##   type          tot prixmed prixmoy surfmed surfmoy prixm2med prixm2moy  part
##   <chr>       <int>   <dbl>   <dbl>   <dbl>   <dbl>     <dbl>     <dbl> <dbl>
## 1 Appartement  9377  145000  174955      45      50      3562      3625  36.5
## 2 Maison      16279  247000  279364      95     102      2652      2806  63.5

1.3 Exploration des prix et des volumes à l’échelon des Iris

1.3.1 Importation des couches géographiques depuis happign

cog <- happign::get_wfs(
  x = NULL,
  apikey = "administratif",
  layer = "ADMINEXPRESS-COG-CARTO.LATEST:commune",
  ecql_filter = "insee_dep = '17'"
)

## Features downloaded : 463

epci_lr <- cog %>%
  filter(siren_epci == "241700434")

# Récupération des Iris par la liste des communes de l'EPCI
liste_com_epci <- list(unique(epci_lr$insee_com))
liste_com_epci <- stringi::stri_join_list(liste_com_epci, sep="','")

IRISLR <- happign::get_wfs(
  x = NULL,
  apikey = "cartovecto",
  layer = "STATISTICALUNITS.IRIS:contours_iris",
  ecql_filter = paste0('insee_com in (\'', liste_com_epci, '\')')
)

## Features downloaded : 67

# Récupération des communes de l'EPCI et des communes attenantes
cog$inter <- st_intersects(cog,st_union(epci_lr), sparse=F)
com_LR_CA <- cog %>%  filter(inter == T) 

## Warning: Using one column matrices in `filter()` was deprecated in dplyr 1.1.0.
## ℹ Please use one dimensional logical vectors instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

liste_com_lrca <- list(unique(com_LR_CA$insee_com))
liste_com_lrca <- stringi::stri_join_list(liste_com_lrca, sep="','")

IRISLR_CA <- happign::get_wfs(
  x = NULL,
  apikey = "cartovecto",
  layer = "STATISTICALUNITS.IRIS:contours_iris",
  ecql_filter = paste0('insee_com in (\'', liste_com_lrca, '\')')
)

## Features downloaded : 80

cog <- cog %>% st_transform(2154)
com_LR_CA <- com_LR_CA %>% st_transform(2154)
IRISLR <- IRISLR %>% st_transform(2154)
IRISLR_CA <- IRISLR_CA %>% st_transform(2154)

1.3.2 Calculer et cartographier le prix moyen au m² par Iris

# Calcul du prix moyen au m² par IRIS dans La Rochelle Agglomération
DVF <- st_as_sf(DVF, coords = c('longitude', 'latitude'), crs = 4326) 
DVF <- st_transform(DVF, 2154)
IRISDVF <- aggregate(dplyr::select(DVF, prix_m2), by = IRISLR, FUN = mean)
IRISDVF <- IRISDVF %>% drop_na()


# Sélection des communes dont le nom sera affiché sur la carte (top 6 selon le nb d'habitants)
liste_noms_LR <- c("La Rochelle", "Aytré", "Puilboreau", "Châtelaillon-Plage", "Lagord", "Périgny")
Selection_Communes_LR <- cog[which(cog$nom %in% liste_noms_LR),] 


# Affichage de la carte
par(oma = c(1,0,0,0))
# Génération du fond de carte
fd1 <- get_tiles(IRISDVF,"CartoDB.PositronNoLabels", 12,crop=T)
fd1 <- project(fd1, "epsg:2154")
plot.new()
mf_raster(fd1)

# Génération de la carte choroplète
mf_map(
  type = "choro",
  x = IRISDVF,
  var = "prix_m2",
  method = "quantiles",
  nbreaks = 5,
  pal = "YlOrRd",
  border = "black",
  leg_title = "Prix moyen \nau m² (Euros)",
  leg_pos = c(399000,6564450),
  add = TRUE)

## The following argument is not relevant when using type = 'choro': method.

# Ajouter les étiquettes des communes selectionnées plus haut, et les autres éléments d'habillage
mf_label(
  Selection_Communes_LR,
  var = "nom",
  halo = TRUE,
  bg = "white",
  r = 0.1,
  cex = 0.8,
  overlap = FALSE)

mf_title("Prix moyen au m² de l'immobilier par Iris sur La Rochelle agglomération et communes attenantes (2014-2022)", cex = 0.8)
mf_credits("IGN et DGFip")
mf_scale(size = 5) 
mf_arrow("topleft")

1.3.3 Calculer et cartographier le nombre de mutations par Iris

# Calcul du nombre de transactions par IRIS
IRISDVF$cpt_vente <- lengths(st_intersects(IRISDVF,DVF))
IRISDVF <- IRISDVF %>% drop_na()

# Affichage de la carte
par(oma = c(1,0,0,0))

# Génération du fond de carte
fd1 <- get_tiles(IRISDVF,"CartoDB.PositronNoLabels", 12,crop=T)
fd1 <- project(fd1, "epsg:2154")
plot.new()
mf_raster(fd1)
plot(st_geometry(IRISDVF), add=TRUE)

# Génération de la carte en symboles proportionnels
mf_map(
  type = "prop",
  x = IRISDVF,
  var = "cpt_vente",
  inches = 0.15,
  col = rgb(0,0.2,1,alpha = 0.5),
  border = rgb(0,0,0,alpha = 0.5),
  leg_title = "Nombre de transactions \npar iris",
  add = TRUE)

# Ajouter les étiquettes des communes selectionnées plus haut, et les autres éléments d'habillage
mf_label(
  Selection_Communes_LR,
  var = "nom",
  halo = TRUE,
  bg = "white",
  r = 0.1,
  cex = 0.8,
  overlap = FALSE)

mf_title("Nombre de mutations par Iris (2014-2022) à La Rochelle agglomération et communes attenantes", cex = 0.8)
mf_credits("IGN et DGFip")
mf_scale(size = 5) 
mf_arrow("topleft")

1.4 Exploration des prix et des volumes à l’échelon des sections cadastrales

1.4.1 Importer la couche des sections cadastrales (DGFip) depuis happign

# Récupération via une liste des communes de l'EPCI 
liste_com <- unique(IRISLR$insee_com)
liste_com <- gsub("^(\\d{2})(.*)", "\\2", liste_com, perl=T)
st <- list(liste_com)
str_com <- stringi::stri_join_list(st, sep="','")
sec_cad <- happign::get_wfs(
  x = NULL,
  apikey = "parcellaire",
  layer = "BDPARCELLAIRE-VECTEUR_WLD_BDD_WGS84G:divcad",
  ecql_filter =  paste0('code_com in (\'', str_com, "') and code_dep = '17'")
)

## Features downloaded : 564

# Même opération en rajoutant les communes attenantes à l'EPCI
liste_com_ca <- unique(IRISLR_CA$insee_com)
liste_com_ca <- gsub("^(\\d{2})(.*)", "\\2", liste_com_ca, perl=T)
st <- list(liste_com_ca)
str_com <- stringi::stri_join_list(st, sep="','")
sec_cad_ca <- happign::get_wfs(
  x = NULL,
  apikey = "parcellaire",
  layer = "BDPARCELLAIRE-VECTEUR_WLD_BDD_WGS84G:divcad",
  ecql_filter =  paste0('code_com in (\'', str_com, "') and code_dep = '17'")
)

## Features downloaded : 848

plot(sec_cad["code_com"])

sec_cad <- st_transform(sec_cad,2154)
sec_cad_ca <- st_transform(sec_cad_ca,2154)

1.4.2 Calculer et cartographier le prix moyen au m² par section cadastrale

# Calcul du prix moyen au m² par section
sec_cadDVF <- sec_cad %>%
  st_join(DVF) %>%
  group_by(id) %>%
  summarise(prix_m2 = round(mean(prix_m2))) %>% 
  na.omit() 

# Affichage de la carte
par(mar = c(1, 0, 0, 0))
# Génération du fond de carte
fd1 <- get_tiles(IRISDVF,"CartoDB.PositronNoLabels", 12,crop=T)
fd1 <- project(fd1, "epsg:2154")
plot.new()
mf_raster(fd1)
plot(st_geometry(sec_cad), add=TRUE)

# Génération de la carte choroplète
mf_map(
  type = "choro",
  x = sec_cadDVF,
  var = "prix_m2",
  method = "quantiles",
  nbreaks = 6,
  pal = "YlOrRd",
  border = rgb(0,0,0,alpha=0.5),
  leg_title = "Prix moyen \nau m² (Euros)",
  add = TRUE)

## The following argument is not relevant when using type = 'choro': method.

# Ajout des étiquettes et autres éléments d'habillage
mf_label(
  Selection_Communes_LR,
  var = "nom",
  halo = TRUE,
  bg = "white",
  r = 0.1,
  cex = 0.8,
  overlap = FALSE)

mf_title("Prix moyen au m² de l'immobilier par section cadastrale à La Rochelle agglo (2014-2022)", cex = 0.8)
mf_credits("IGN et DGFip")
mf_scale(size = 5) 
mf_arrow("topleft")

1.4.3 Calculer et cartographier le nombre de mutations par section cadastrale

# Calcul du nombre de mutations par section
sec_cadDVF <- sec_cadDVF %>%
  mutate(Nbmutations = lengths(st_intersects(sec_cadDVF, DVF)))

# Affichage de la carte
par(mar = c(1, 0, 0, 0))
# Génération du fond de carte
fd1 <- get_tiles(IRISDVF,"CartoDB.PositronNoLabels", 12,crop=T)
fd1 <- project(fd1, "epsg:2154")
mf_raster(fd1)
plot(st_geometry(sec_cad), border= rgb(0,0,0,alpha=0.3), add=T)

# Génération de la carte en symboles proportionnels
propSymbolsLayer(x = sec_cadDVF, 
                 var = "Nbmutations", 
                 col = rgb(0,0.2,1,alpha = 0.5),
                 border = "white",  
                 lwd = 0.01,
                 inches = 0.1, 
                 add = TRUE,
                 legend.pos = "n") 

# Ajout des étiquettes et autres éléments d'habillage
layoutLayer(title = "Nombre de mutations par section cadastrale à La Rochelle Agglomération (2014-2022)", source = "         IGN et DGFip", north = TRUE, horiz = FALSE, tabtitle = TRUE, frame = FALSE, col = "black", coltitle = "white")

legendCirclesSymbols(
  pos = "bottomleft",
  title.txt = "Nombre de mutations DVF",
  title.cex = 0.8,
  cex = 1,
  border = "white",
  lwd = 0.01,
  values.cex = 0.6,
  var = c(min(sec_cadDVF$Nbmutations), max(sec_cadDVF$Nbmutations)),
  inches = 0.08,
  col = rgb(0,0.2,1,alpha = 0.5),
  frame = FALSE,
  values.rnd = 0,
  style = "e"
)

# Ajouter les étiquettes des communes selectionnées plus haut
labelLayer(Selection_Communes_LR, txt = "nom", halo = TRUE, bg = "white", r = 0.1, cex = 0.8, pos = 3, font = 3, offset = 0, overlap = TRUE)

1.5 Lissage spatial (échelle métropolitaine)

1.5.1 Lissage spatial des prix de l’immobilier au m² calculés à partir des valeurs agrégées par Iris

1.5.1.1 Calcul des distances au plus proche voisin et de l’auto-corrélation spatiale des prix au m²

# pour calculer les pseudo-centroides des Iris 
IRISLR.Centroids <- st_point_on_surface(IRISLR)

# Identification du plus proche voisin de chaque Iris
listPPV_iris_cen <- knearneigh(IRISLR.Centroids, k = 1) 

# Conversion de l'objet knn en objet nb
PPV_iris <- knn2nb(listPPV_iris_cen, row.names = IRISLR.Centroids$code_iris) 

distPPV_iris <- nbdists(PPV_iris, IRISLR.Centroids) # distance entre plus proches voisins
print(as.data.frame(t(as.matrix(summary(unlist(distPPV_iris)))))) # résumé statistique

##       Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
## 1 331.1623 543.4473 1081.793 1339.615 1882.393 4854.631

# des distances entre IRIS

hist(unlist(distPPV_iris), breaks = 20,
     main = "Distance au plus proche voisin",
     col = "black", border = "white", xlab = "Distance", ylab = "Fréquence")

# pour convertir les Iris en objet nb
IRIS.nb <- poly2nb(pl = IRISDVF,
                   snap = 50,
                   queen = TRUE) # voisinage au sens large sommet ou frontière communes

#calcul du test de Moran
moran.test(IRISDVF$prix_m2, listw = nb2listw(IRIS.nb))

## 
##  Moran I test under randomisation
## 
## data:  IRISDVF$prix_m2  
## weights: nb2listw(IRIS.nb)    
## 
## Moran I statistic standard deviate = 8.1923, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.620170543      -0.015151515       0.006014213

Relation supérieur (Moran I) à 0 doc autocorrélation positive. Z-score = 8.1923, en dehors de l’intervalle [-1,96 ; 1,96] pour \(\alpha\) = 0.95~%. On lisse quand on a de l’autocorrélation spatiale, sinon c’est déconseillé. Ici, l’indice de Moran nous autorise à représenter les données par un lissage spatial

1.5.2 Carte lissée des prix de l’immobilier au m² calculés à partir des valeurs agrégées par Iris sur l’ensemble de La Rochelle Agglomération (et communes attenantes)

Nous décidons de réaliser le lissage en prenant en compte les communes attenantes afin de limiter l’effet de bord.

# pour définir le contour de La Rochelle agglomération (et communes attenantes) comme emprise pour le lissage
EmpriseLR <- as.owin(st_geometry(com_LR_CA))
IRISDVF_CA <- aggregate(dplyr::select(DVF, prix_m2), by = IRISLR_CA, FUN = mean)
IRISDVF_CA <- IRISDVF_CA %>% drop_na()

# pour récupérer les coordonnées du centroïde des IRIS
pts <- st_coordinates(st_point_on_surface(st_geometry(IRISLR_CA)))

# pour créer un objet ppp (format spatstat) et y intégrer dedans l'emprise et les valeurs à lisser (prix au m²)
IRISDVF_CA.ppp <- ppp(pts[,1], pts[,2], window = EmpriseLR, marks = IRISDVF_CA$prix_m2)

# pour calculer la surface lissée (rayon lissage : 1 km et resolution spatiale de l'image : 1 ha)
cartelissee <- Smooth(IRISDVF_CA.ppp, sigma = 1000, weights = IRISDVF_CA.ppp$marks, eps = 100)

# Conversion de la surface lissée au format raster et découpage selon l'emprise
# de l'EPCI
# Conversion de la surface lissée au format raster
cartelissee.raster <- terra::rast(cartelissee)
crs(cartelissee.raster)  <- "epsg:2154"
# Convertir epci_lr en spatVector 
mask <- vect(st_transform(epci_lr,2154))
masked <- terra::mask(cartelissee.raster, mask)

# Paramétrage des marges de la fenêtre pour maximiser l'emprise de la carte
par(mar = c(1, 0, 0, 0))

# Pour afficher la surface lissée et définir l'habillage
# Calcul des seuils
bks <- unique(getBreaks(values(masked), method = "q6"))

# Création d'une palette de couleurs (double gradation harmonique)
cols <- c("#2a9c4e", "#77c35c", "#c4e687", "#ffffc0", "#fec981", "#dc292c")

fd1 <- get_tiles(IRISLR_CA,"CartoDB.PositronNoLabels", 12, crop=T)
fd1 <- project(fd1, "epsg:2154")
mf_raster(fd1)

# Affichage de la carte lissée et du contour des Iris
plot(masked, breaks = bks, col = cols, add = T, legend = F)
plot(epci_lr$geometry, border = rgb(1,1,1, alpha = 0.5), lwd = 0.5, add = T)
plot(IRISDVF$geometry, border = "grey60", lwd = 0.05, lty = 3, add = T)

legendChoro(
  pos = "bottomleft",
  title.txt = "Prix moyen \nau m² (Euros)",
  breaks = bks, 
  nodata = FALSE,
  values.rnd = -1,
  border = "white",
  col = cols
)

# Habillage de la carte
layoutLayer(title = "Prix moyen au m² de l'immobilier à La Rochelle Agglo (2014-2022)", 
            author = "          Sources : IGN et DGFip - Rayon : 1000 m,\n          résolution : 100m - A partir des Iris", 
            scale = 5, frame = TRUE, col = "black", coltitle = "white", 
            north = TRUE, tabtitle = TRUE, horiz = FALSE)

# Ajouter les étiquettes des communes sélectionnées plus haut
labelLayer(Selection_Communes_LR, txt = "nom", halo = TRUE, bg = "white", r = 0.05, cex = 0.6, pos = 3, font = 3, offset = 0, overlap=F)

1.5.3 Lissage spatial du prix de l’immobilier au m² calculés à partir des valeurs agrégées par section cadastrale

1.5.3.1 Calcul des distances au plus proche voisin et de l’auto-corrélation spatiale des prix au m²

# Récupération des valeurs de prix_m2 par section cadastrale
sec_cad_caDVF <- sec_cad_ca %>%
  st_join(DVF) %>%
  group_by(id) %>%
  summarise(prix_m2 = round(mean(prix_m2)))
sec_cad_caDVF <- sec_cad_caDVF %>% na.omit()

# pour calculer les pseudo-centroides des sections cadastrales (indispensable pour le calcul sur les plus proches voisins)
sec_cad.Centroids <- st_point_on_surface(sec_cad_caDVF)

# Calcul sur les plus proches voisins
listPPV <- knearneigh(sec_cad.Centroids, k = 1) # pour connaître le plus proche voisin de chaque section cadastrale
PPV <- knn2nb(listPPV, row.names = sec_cad.Centroids$id) # pour convertir l'objet knn en objet nb
distPPV <- nbdists(PPV, sec_cad.Centroids) # pour connaître la distance entre plus proches voisins

print(as.data.frame(t(as.matrix(summary(unlist(distPPV))))))

##       Min. 1st Qu.  Median     Mean  3rd Qu.     Max.
## 1 126.1835 405.637 514.648 581.8634 712.9702 3162.455

hist(unlist(distPPV), breaks = 20,
     main = "Distance au plus proche voisin",
     col = "black", border = "white", xlab = "Distance", ylab = "Fréquence")

# pour convertir les sections en objet nb
sec_cad_ca.nb <- poly2nb(pl = sec_cad_caDVF,
                 snap = 980,
                 queen = TRUE) #Nous augmentons le snap à 980, car l'indice 
#de Moran ne peut pas être calculé avec une valeur inférieure

#calcul du test de Moran
moran.test(sec_cad_caDVF$prix_m2, listw = nb2listw(sec_cad_ca.nb))

## 
##  Moran I test under randomisation
## 
## data:  sec_cad_caDVF$prix_m2  
## weights: nb2listw(sec_cad_ca.nb)    
## 
## Moran I statistic standard deviate = 35.789, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.6004203228     -0.0016977929      0.0002830445

1.5.3.2 Carte lissée des prix de l’immobilier au m² calculés à partir des valeurs agrégées par section cadastrale sur l’ensemble de Rennes Métropole

# pour récupérer les coordonnées du centroïde des sections
pts <- st_coordinates(st_point_on_surface(st_geometry(sec_cad_caDVF)))

# pour créer un objet ppp (format spatstat) et y intégrer dedans l'emprise et les valeurs à lisser (prix au m²)
sec_cad.ppp <- ppp(pts[,1], pts[,2], window = EmpriseLR, marks = sec_cad_caDVF$prix_m2)

## Warning: 4 points were rejected as lying outside the specified window

# pour calculer la surface lissée (rayon lissage : 1 km et resolution spatiale de l'image : 1 ha)
cartelissee <- Smooth(sec_cad.ppp, sigma = 1000, weights = sec_cad.ppp$marks, eps = 100)

# Conversion de la surface lissée au format raster
cartelissee.raster <- terra::rast(cartelissee)
crs(cartelissee.raster)  <- "epsg:2154"

# Convertir epci_lr en spatVector 
mask <- vect(st_transform(epci_lr,2154))
masked <- terra::mask(cartelissee.raster, mask)

# Paramétrage des marges de la fenêtre pour maximiser l'emprise de la carte
par(mar = c(1, 0, 0, 0))

# pour afficher la surface lissée et définir l'habillage
# Calcul des seuils
bks <- unique(getBreaks(values(masked), method = "q6"))

# Création d'une palette de couleurs (double gradation harmonique)
cols <- c("#2a9c4e", "#77c35c", "#c4e687", "#ffffc0", "#fec981", "#dc292c")

fd1 <- get_tiles(sec_cad_caDVF,"CartoDB.PositronNoLabels", 12,crop=T)
fd1 <- project(fd1, "epsg:2154")

mf_raster(fd1)

# Affichage de la carte lissée et du contour des Iris
plot(masked, breaks = bks, col = cols, add = T, legend = F)
plot(sec_cad_ca$geometry, border = rgb(0,0,0,alpha = 0.5), lwd = 0.03, lty = 3, add = T)
plot(com_LR_CA$geometry, border = "white", lwd = 0.5, add = T)

legendChoro(
  pos = "bottomleft",
  title.txt = "Prix moyen \nau m² (Euros)",
  breaks = bks, 
  nodata = FALSE,
  values.rnd = -1,
  border = "white",
  col = cols
)

# Habillage de la carte
layoutLayer(title = "Prix moyen au m² de l'immobilier à La Rochelle Agglo (2014-2022)", 
            author = "          Sources : IGN et DGFip - Rayon : 1000 m,\n          résolution : 100m - A partir des Iris", 
            scale = 5, frame = TRUE, col = "black", coltitle = "white", 
            north = TRUE, tabtitle = TRUE, horiz = FALSE)

# Ajouter les étiquettes des communes sélectionnées plus haut
labelLayer(Selection_Communes_LR, txt = "nom", halo = TRUE, bg = "white", r = 0.05, cex = 0.6, pos = 3, font = 3, offset = 0, overlap=F)

1.5.4 Lissage spatial des prix de l’immobilier au m² calculés directement à partir des mutations immobilières sur l’ensemble de La Rochelle agglomération (et communes attenantes)

1.5.4.1 Calcul de la distance au plus proche voisin et de l’auto-corrélation spatiale des prix au m²

# Calcul sur les plus proches voisins
listPPV <- knearneigh(DVF, k = 1) #pour connaître le plus proche voisin de chaque mutation

PPV <- knn2nb(listPPV, row.names = DVF$id_vente) # pour convertir l'objet knn en objet nb

distPPV <- nbdists(PPV, DVF) # pour connaître la distance entre plus proches voisins

print(as.data.frame(t(as.matrix(summary(unlist(distPPV))))))

##   Min. 1st Qu.   Median     Mean 3rd Qu.     Max.
## 1    0       0 6.784495 17.69769 26.6294 2459.984

hist(unlist(distPPV), nclass = 500,
     main = "Distance au plus proche voisin",
     col = "black", border = "white", xlab = "Distance", ylab = "Fréquence", xlim = c(0,150))

#calcul du test de Moran
moran.test(DVF$prix_m2, listw = nb2listw(PPV))

## 
##  Moran I test under randomisation
## 
## data:  DVF$prix_m2  
## weights: nb2listw(PPV)    
## 
## Moran I statistic standard deviate = 61.473, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      4.690504e-01     -3.897876e-05      5.823001e-05

1.5.4.2 Carte lissée des prix de l’immobilier au m² calculés directement à partir des mutations immobilières sur l’ensemble de La Rochelle agglomération (et communes attenantes)

# Extraction des coordonnées spatiales
coordinates <- st_coordinates(DVF)

# pour créer un objet ppp (format spatstat) et y intégrer dedans l'emprise et les valeurs à lisser (prix moyen au m²)
DVF.ppp <- ppp(x = coordinates[, 1], y = coordinates[, 2], window = EmpriseLR, marks = DVF$prix_m2)

# pour calculer la surface lissée (rayon lissage : 0,5 km et résolution spatiale de l'image : 1 ha) --> calcul long
cartelissee <- Smooth(DVF.ppp, sigma = 500, weights = DVF.ppp$marks, eps = 100)

# Conversion de la surface lissée au format raster
cartelissee.raster <- terra::rast(cartelissee)
crs(cartelissee.raster) <- st_crs(DVF)$srid

# découpage de la surface lissée sur l'emprise Rennes Métropole
cartelissee.raster <- mask(cartelissee.raster, st_transform(epci_lr,2154)) 

# Paramétrage des marges de la fenêtre pour maximiser l'emprise de la carte
par(mar = c(1, 0, 0, 0))

# pour afficher la surface lissée et définir l'habillage de la carte
# Calcul des seuils
bks <- unique(getBreaks(values(cartelissee.raster), method = "q6"))

# Création d'une palette de couleurs (double gradation harmonique)
cols <- c("#2a9c4e", "#77c35c", "#c4e687", "#ffffc0", "#fec981", "#dc292c")

fd1 <- get_tiles(DVF,"CartoDB.PositronNoLabels", 12,crop=T)
fd1 <- project(fd1, "epsg:2154")

mf_raster(fd1)

# Affichage de la carte lissée et du contour des Iris
plot(cartelissee.raster, breaks = bks, col = cols, add = T, legend = F)
plot(com_LR_CA$geometry, border = "white", lwd = 0.5, add = T)

legendChoro(
  pos = "bottomleft",
  title.txt = "Prix moyen \nau m² (Euros)",
  breaks = bks, 
  nodata = FALSE,
  values.rnd = -1,
  border = "white",
  col = cols
)

# Habillage de la carte
layoutLayer(title = "Prix moyen au m² de l'immobilier à La Rochelle Agglo (2014-2022)", 
            author = "          Sources : IGN et DGFip - Rayon : 500 m,\n          résolution : 100m - A partir des Iris", 
            scale = 5, frame = TRUE, col = "black", coltitle = "white", 
            north = TRUE, tabtitle = TRUE, horiz = FALSE)

# Ajouter les étiquettes des communes sélectionnées plus haut
labelLayer(Selection_Communes_LR, txt = "nom", halo = TRUE, bg = "white", r = 0.05, cex = 0.6, pos = 3, font = 3, offset = 0, overlap=F)

1.6 Elaboration d’une typologie

Une typologie est établie à partir d’indicateurs immobiliers et une cartographie des sous-marchés associés est réalisée.

1.6.1 Transfert d’attributs des Iris vers les mutations (jointure spatiale)

MutationsIRIS <- DVF %>%
  st_join(IRISLR) %>%
  drop_na()

MutationsIRIS <- as.data.frame(MutationsIRIS)

1.6.2 Création d’un tableau avec les variables par Iris à soumettre à la CAH

MutationsIRIS <- MutationsIRIS %>%
  mutate(prixm2 = prix / surface)

IRISDVFClassif1 <- MutationsIRIS %>%
  group_by(code_iris) %>%
  summarise(Nbtransactions = n(),
            code_iris = first(code_iris),
            Prixmoyen = round(mean(prix)),
            Prixm2moyen = round(mean(prixm2)),
            Surfacemoyenne = round(mean(surface)),
            PropMaison = length(type[type == "Maison"]) / Nbtransactions * 100,
            PropAppart = length(type[type == "Appartement"]) / Nbtransactions * 100)

IRISDVFClassif <- data.frame(IRISDVFClassif1[, c("Nbtransactions", "Prixmoyen", "Prixm2moyen", "Surfacemoyenne", "PropMaison", "PropAppart")])

Les variables sont centrées réduites afin de standardiser les données.

IRISDVFClassifScale <- scale(IRISDVFClassif)

1.6.3 Mise en oeuvre de la CAH

CAHIRIS <- agnes(
   IRISDVFClassifScale,
   metric = "euclidean",
   method = "ward"
)

1.6.4 Graphiques des gains d’inertie inter-classe

sortedHeight <- sort(CAHIRIS$height, decreasing = TRUE)
relHeight <- sortedHeight/ sum(sortedHeight) * 100
cumHeight <- cumsum(relHeight)

library(Factoshiny)

## Le chargement a nécessité le package : FactoMineR

## Le chargement a nécessité le package : shiny

## Le chargement a nécessité le package : FactoInvestigate

library(FactoMineR)
pca <- PCA(IRISDVFClassifScale, ncp = Inf, scale.unit = FALSE, graph = FALSE)
cah <- HCPC(pca, graph = FALSE)
plot(cah, choice = "bar")

# Gain d'inertie
barplot(
  relHeight[1:30],
  names.arg = seq(1, 30, 1),
  col = "black", border = "white",
  xlab = "Noeuds", ylab = "Part de l'inertie totale (%)"
)

barplot(
  cumHeight[1:30],
  names.arg = seq(1, 30, 1),
  col = "black", border = "white",
  xlab = "Nombre de classes", ylab = "Part de l'inertie totale (%)"
)

1.6.5 Arbre de classification hiérarchique

dendroCSP <- as.dendrogram(CAHIRIS)
plot(dendroCSP, leaflab = "none")

1.6.6 Partition (en n classes)

clusIRIS <- cutree(CAHIRIS, k = 3)
IRISCluster <- IRISDVFClassif1
IRISCluster$CLUSIMMO <- factor(
   clusIRIS,
   levels = 1:3,
   labels = paste("Classe", 1:3)
)

1.6.7 Création d’un tableau avec les caractéristiques des groupes

RecapCAHIRIS <- IRISCluster %>%
  group_by(CLUSIMMO) %>%
  summarise(
    NB= n(),
    NbTransac = round(mean(Nbtransactions)),
    Prixmoyen = round(mean(Prixmoyen)),
    Prixm2 = round(mean(Prixm2moyen)),
    Surface = round(mean(Surfacemoyenne)),
    PropMaison = mean(PropMaison),
    PropAppart = mean(PropAppart)
  )

print(RecapCAHIRIS)

## # A tibble: 3 × 8
##   CLUSIMMO    NB NbTransac Prixmoyen Prixm2 Surface PropMaison PropAppart
##   <fct>    <int>     <dbl>     <dbl>  <dbl>   <dbl>      <dbl>      <dbl>
## 1 Classe 1    37       293    271011   2722     102       90.1       9.92
## 2 Classe 2    18       182    197212   2611      77       66.9      33.1 
## 3 Classe 3    12       713    217482   3775      60       23.7      76.3

1.6.8 Graphique des écarts à la moyenne

SyntheseCAHIRIS <- RecapCAHIRIS %>% mutate(
  nbtransacmoy = mean(IRISDVFClassif$Nbtransactions),
  surfacemoy = mean(IRISDVFClassif$Surfacemoyenne),
  prixmoy = mean(IRISDVFClassif$Prixmoyen),
  prixm2moyen = mean(IRISDVFClassif$Prixm2moyen),
  propmaisonmoyen = mean(IRISDVFClassif$PropMaison),
  propappartmoyen = mean(IRISDVFClassif$PropAppart),
  NbMutations = (NbTransac - nbtransacmoy) / nbtransacmoy * 100,
  Prix = (Prixmoyen - prixmoy) / prixmoy * 100,
  Prixm2 = (Prixm2 - prixm2moyen) / prixm2moyen * 100,
  Surface = (Surface - surfacemoy) / surfacemoy * 100,
  PropMaison = (PropMaison - propmaisonmoyen) / propmaisonmoyen * 100,
  PropAppart = (PropAppart - propappartmoyen) / propappartmoyen * 100
)

# pour créer un nouveau tableau avec une sélection de colonnes à faire apparaître sur le graphique
SyntheseCAHIRIS <- data.frame(
  SyntheseCAHIRIS[,
    c("CLUSIMMO", "NbMutations", "Surface", "Prix", "Prixm2", "PropMaison", "PropAppart")
  ]
)

gather <- SyntheseCAHIRIS %>%
  gather(key = variable, value = "value", NbMutations:PropAppart)

ggplot(gather, aes(x = variable, y = value, fill = CLUSIMMO)) +
  geom_bar(stat = "identity") +
  coord_flip() +
  scale_fill_manual(values = c("#2A9D8F","#264653","#fcbf49","#e45c3a","#0074D9")) +
  ylab("Variation par rapport à la moyenne métropolitaine (%)") +
  theme_bw() +
  theme(legend.position = "none") +
  facet_wrap(~CLUSIMMO, ncol = 1)

1.6.9 Carte de la typologie montrant les sous-marchés immobiliers métropolitains

# pour joindre le résultat de la typologie dans la couche des IRIS
IRISDVFCAH <- left_join(IRISLR, IRISCluster, by = "code_iris")

par(mar = c(0, 0, 1.2, 2))

# Affichage des communes de RM et alentours en arrière-plan et centrage de la carte sur Rennes Métropole
plot(
   st_geometry(IRISLR),
   lwd = 0.1,
   col = "grey95",
   border = "white",
   bg = "#B5D0D0"
   # xlim = c(st_bbox(IRISRM)[1], st_bbox(IRISRM)[3]), ylim = c(st_bbox(IRISRM)[2], st_bbox(IRISRM)[4])
)

# Affichage de la typologie et du contour des IRIS
unique(IRISDVFCAH$CLUSIMMO)

## [1] Classe 1 Classe 2 Classe 3
## Levels: Classe 1 Classe 2 Classe 3

IRISDVFCAH

## Simple feature collection with 67 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 372869.8 ymin: 6552942 xmax: 398644.9 ymax: 6581807
## Projected CRS: RGF93 v1 / Lambert-93
## # A tibble: 67 × 17
##    id            insee_com nom_com iris  code_iris nom_iris typ_iris layer path 
##    <chr>         <chr>     <chr>   <chr> <chr>     <chr>    <chr>    <chr> <chr>
##  1 IRIS____0000… 17291     Puilbo… 0000  172910000 Puilbor… Z        CONT… R:/c…
##  2 IRIS____0000… 17300     La Roc… 0505  173000505 Mireuil… H        CONT… R:/c…
##  3 IRIS____0000… 17300     La Roc… 0202  173000202 La Gene… H        CONT… R:/c…
##  4 IRIS____0000… 17274     Périgny 0103  172740103 Périgny… H        CONT… R:/c…
##  5 IRIS____0000… 17420     Salles… 0000  174200000 Salles-… Z        CONT… R:/c…
##  6 IRIS____0000… 17300     La Roc… 0402  173000402 La Pall… H        CONT… R:/c…
##  7 IRIS____0000… 17153     Esnand… 0000  171530000 Esnandes Z        CONT… R:/c…
##  8 IRIS____0000… 17300     La Roc… 1001  173001001 Les Min… H        CONT… R:/c…
##  9 IRIS____0000… 17200     Lagord  0103  172000103 Ouest    H        CONT… R:/c…
## 10 IRIS____0000… 17300     La Roc… 0903  173000903 Villene… H        CONT… R:/c…
## # ℹ 57 more rows
## # ℹ 8 more variables: geometry <MULTIPOLYGON [m]>, Nbtransactions <int>,
## #   Prixmoyen <dbl>, Prixm2moyen <dbl>, Surfacemoyenne <dbl>, PropMaison <dbl>,
## #   PropAppart <dbl>, CLUSIMMO <fct>

typoLayer(
  x = IRISDVFCAH,
  var="CLUSIMMO",
  col = c("#2A9D8F","#264653","#fcbf49","#e45c3a","#0074D9"),
  lwd = 0.05,
  border = "grey70",
  legend.pos = "bottomleft",
  legend.title.txt = "Sous-marchés \nimmobiliers",
  legend.nodata = "Aucune mutation",
  add = TRUE)

layoutLayer(title = "Sous-marchés immobiliers dans Rennes Métropole à l'échelon des IRIS (2014-2019)",
            author = "          Sources : IGN et DGFip - Typologie obtenue par CAH",
            scale = 5, frame = TRUE, col = "black", coltitle = "white",
            north = TRUE, tabtitle = TRUE, horiz = FALSE)

# Ajouter les étiquettes des communes sélectionnées plus haut
# FIXME: ça marche pas
labelLayer(Selection_Communes_LR, txt = "NOM_COM", halo = TRUE, bg = "white", r = 0.05, cex = 0.6, pos = 3, font = 3, offset = 0)

1.6.10 Export des Iris contenant la typologie issue de la CAH au format geopackage

Export de classif.

st_write(IRISDVFCAH, "data/IRISDVFCAH.gpkg", append = FALSE)

## Deleting layer `IRISDVFCAH' using driver `GPKG'
## Writing layer `IRISDVFCAH' to data source `data/IRISDVFCAH.gpkg' using driver `GPKG'
## Writing 67 features with 16 fields and geometry type Multi Polygon.

1.7 Lissage spatial des prix au m² des appartements entre 2014 et 2022

Zoom sur la commune de La Rochelle

1.7.1 Préparation du jeu de données

# pour ne garder que la commmune de La Rochelle
cog_roch <- cog %>%
  filter(nom == "La Rochelle")

# pour ne garder que les mutations correspondant
# aux appartements dans La Rochelle Agglomération (dans un premier temps)
dvf_appt_roch <- DVF %>%
  filter(codecommune == 17300 & type == "Appartement")

# pour créer une liste d'années, dans l'ordre
dvf_appt_roch <- dvf_appt_roch[order(dvf_appt_roch$annee, decreasing = FALSE), ]
ListAnnees <- unique(dvf_appt_roch$annee)

# pour définir comme emprise pour le lissage, l'étendue
# de La Rochelle avec une zone tampon de 0,5 km autour afin
# d'englober les mutations limitrophes et réduire l'effet de bord
bbox_roch <- as.owin(st_buffer(st_geometry(cog_roch), 500))

# Création d'une palette de couleurs
cols <- c(
  "#fff7f3",
  "#fee1d3",
  "#fcbea5",
  "#fc9777",
  "#fb7050",
  "#f14432",
  "#d42020",
  "#ad1016",
  "#67000d",
  "#480000",
  "#300000")

1.7.2 Cartes lissées

Série de cartes lissées montrant l’évolution des prix au m² des appartements calculés directement à partir des mutations de la Rochelle

# Paramétrage des marges pour insérer le titre général et les titres de chaque carte
par(
  oma = c(3.5, 0, 3, 0) + 0.1,
  mar = c(0, 0.5, 1.2, 0.5)
)

# pour découper la fenêtre en 3 lignes et 3 colonnes (9 années)
par(mfrow = c(3, 3))

# boucle pour produire et cartographier une surface lissée par année
for (i in ListAnnees){
  print(i)

  # Recuperation des jeux de données par année
  IRISDVF_CA <- DVF %>%
    filter(annee == i) %>%
    drop_na(prix_m2)

  print(paste("Number of mut =", length(IRISDVF_CA$prix_m2)))
  print(paste("Number of mut =", mean(IRISDVF_CA$prix_m2)))
  print(paste("Number of mut =", sum(IRISDVF_CA$prix_m2)))

  IRISDVF_CA <- aggregate(
    dplyr::select(IRISDVF_CA, prix_m2), by = IRISLR_CA, FUN = mean
  )
  IRISDVF_CA <- IRISDVF_CA[which(!is.na(IRISDVF_CA$prix_m2)), ]

  # pour récupérer les coordonnées du centroïde des IRIS
  pts <- st_coordinates(st_point_on_surface(st_geometry(IRISDVF_CA)))
  # pour créer un objet ppp (format spatstat) et y intégrer dedans
  # l'emprise et les valeurs à lisser (prix au m²)
  IRISDVF_CA.ppp <- ppp(pts[,1], pts[,2], window = EmpriseLR, marks = IRISDVF_CA$prix_m2)

  # pour calculer la surface lissée (rayon lissage : 250m et resolution 
  # spatiale de l'image : 100 m)
  cartelissee <- Smooth(IRISDVF_CA.ppp, sigma = 250, weights = IRISDVF_CA.ppp$marks, eps = 100)

  # Conversion de la surface lissée au format raster
  cartelissee.raster <- rast(cartelissee)
  crs(cartelissee.raster) <- st_crs(cog_roch)$srid # pour spécifier un SCR à l'objet raster
  cartelissee.raster <- mask(cartelissee.raster, cog_roch) # découpage sur emprise La Rochelle

  # Définition manuelle des seuils
  bks <- c(min(values(cartelissee.raster), na.rm = TRUE), 1500, 2000, 2200, 2400, 2800, 3000, 3200, 3500, 4000, max(values(cartelissee.raster), na.rm = TRUE))

  # reclassification de la surface lissée
  cartelissee.reclass <- classify(cartelissee.raster, rcl = bks)

  # vectorisation de la surface reclassée
  cartelissee.vecteur <- st_as_sf(as.polygons(cartelissee.reclass))

  # trier les valeurs
  cartelissee.vecteur <- cartelissee.vecteur[order(cartelissee.vecteur$lyr.1, decreasing = FALSE), ]

  # Tracer la carte
  plot(st_geometry(cog_roch), border = "white", bg = "grey90")

  typoLayer(
    x = cartelissee.vecteur, var = "lyr.1", col = cols, lwd = 0.1,
    border = NA, legend.pos = "n", add = TRUE
  )

  plot(
    st_geometry(IRISLR),
    border = adjustcolor("white", alpha.f = 0.5),
    lty = 3, lwd = 0.01, add = TRUE
  )

  title(main = paste("", i, sep = ""))

}

## [1] 2014
## [1] "Number of mut = 2379"
## [1] "Number of mut = 2822.79118247799"
## [1] "Number of mut = 6715420.22311514"

## [1] 2015
## [1] "Number of mut = 2666"
## [1] "Number of mut = 2764.65239661806"
## [1] "Number of mut = 7370563.28938374"

## [1] 2016
## [1] "Number of mut = 3099"
## [1] "Number of mut = 2763.4410373335"
## [1] "Number of mut = 8563903.77469651"

## [1] 2017
## [1] "Number of mut = 3733"
## [1] "Number of mut = 2900.45902116142"
## [1] "Number of mut = 10827413.5259956"

## [1] 2018
## [1] "Number of mut = 3827"
## [1] "Number of mut = 3005.90376906455"
## [1] "Number of mut = 11503593.72421"

## [1] 2019
## [1] "Number of mut = 4093"
## [1] "Number of mut = 3216.62295294932"
## [1] "Number of mut = 13165637.7464215"

## [1] 2020
## [1] "Number of mut = 2199"
## [1] "Number of mut = 3457.52982466441"
## [1] "Number of mut = 7603108.08443704"

## [1] 2021
## [1] "Number of mut = 2409"
## [1] "Number of mut = 3702.98092072082"
## [1] "Number of mut = 8920481.03801646"

## [1] 2022
## [1] "Number of mut = 1251"
## [1] "Number of mut = 4000.58725714111"
## [1] "Number of mut = 5004734.65868353"

barscale(
  lwd = 1.5, cex = 0.6, pos = "bottomright", style = "pretty", unit = "m"
)
cartography::north(pos = "topleft")

# Pour afficher le titre principal et la source
mtext(
  "Prix moyen en Euros au m² des appartements à La Rochelle de 2014 à 2022",
  cex = 1, side = 3, line = 1, adj = 0.5, outer = TRUE
)
mtext(
  "   Sources : IGN et DGFip - Rayon de lissage : 250 m\n   et résolution : 100 m", 
  side = 1, line = 1, adj = 0, cex = 0.6, outer = TRUE
)

# Affichage de la légende (dans une nouvelle fenêtre graphique)
par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
plot.new()

legendChoro(
  pos = "bottomright",
  title.txt = "",
  breaks = bks,
  nodata = FALSE,
  values.rnd = -1,
  col = cols,
  border = "white",
  horiz = TRUE
)